NX=32;
NY=64;
    Pic_Num=26000;
    v_i=0.05;
    v_o=0.2;
    %v_o=v_i;
    G=1.0e-6;
    a=15.5;
    
   dat1=load(strcat('LBM_v_',int2str(Pic_Num),'.out'));
   dat2=load(strcat('LBM_rho_',int2str(Pic_Num),'.out'));
   dat3=load(strcat('LBM_rhoR_',int2str(Pic_Num),'.out'));
   %dat4=load(strcat('matrix_2_',int2str(Pic_Num),'.dat'));
   %dat3=load(strcat('LBM_velocity_poi_y',int2str(Pic_Num),'.dat'));
   %dat5=load(strcat('LBM_velocity_poi_x',int2str(Pic_Num),'.dat'));
   

   
   u=[];
   %L=NY/2-0.5;
   L=(NY-1)/2;
   x=[-(NY-1)/2-0.5:(NY-1)/2+0.5];
   for i=1:NY+1
       if (abs(x(i))<a )
           u(i)=G/(2*v_o)*(L^2-a^2)+G/(2*v_i)*(a^2-x(i)^2);
       else
           u(i)=G/(2*v_o)*(L^2-x(i)^2);
       end
   %         u(i)=G/(2*v_o)*(L^2-x(i)^2);
   end
   
   figure
   plot(u,'r');
   hold on
   plot (dat1,'bx');
   hold off
 
   figure
   plot (dat2,'r');
   hold on
   plot (dat3,'b');
   hold off